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We study a quantum version of the three-state Potts model that includes as special cases the effective models 
of bosons and fermions on the square lattice in the Mott insulating limit. It can be viewed as a model of quantum 
permutations with amplitudes J\\ and J± for identical and different colors, respectively. For Jy — J± > it is 
equivalent to the SU(3) Heisenberg model, which describes the Mott insulating phase of 3-color fermions, while 
the parameter range J± < min(0, —J\\ ) can be realized in the Mott insulating phase of 3-color bosonic atoms. 
Using linear flavor wave theory, infinite projected entangled-pair states (iPEPS), and continuous-time quantum 
Monte-Carlo simulations, we construct the full T = phase diagram, and we explore the T 7^ properties for 
J±_ < 0. For dominant antiferromagnetic Jy interactions, a three-sublattice long-range ordered stripe state is 
selected out of the ground state manifold of the antiferromagnetic Potts model by quantum fluctuations. Upon 
increasing | J^|, this state is replaced by a uniform superfluid for J± < 0, and by an exotic three-sublattice 
superfluid followed by a two-sublattice superfluid for J± > 0. The transition out of the uniform superfluid (that 
can be realized with bosons) is shown to be a peculiar type of Kosterlitz-Thouless transition with three types of 
elementary vortices. 

PACS numbers: 67.85.-d, 75.10Jm, 02.70.-c 



I. INTRODUCTION 

The possibility to realize N-color Fermi- or Bose-Hubbard 
models with ultra-cold atoms in optical lattices has attracted 
increasing interest recently. ~ With fermionic alkaline-earth 
atoms (and with the alkaline-earth like atomYb), it is possible 
to realize systems with up to N=10 different flavors (or colors) 
of fermions by exploiting different nuclear spin states, while 
3-color bosonic systems can be created with bosonic alkali 
atoms with a hyperfine spin F = 1, e.g. ^^Na, ^^K, and ^^Rb. 
A good starting point to describe these systems is the SU(N) 
Hubbard model given by 

{ij),a 



- it vanishes if the scattering lengths of the spin-0 and spin-2 
channels are equal - and we will not consider it in the present 
paper. 

At filling 1/A^ (one particle per site) and for large enough 
repulsion the system is in a Mott insulating phase, and color 
fluctuations can be described by the effective Hamiltonian: 
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where (ij) denotes a sum over all pairs of nearest neighbors, 
and the operators cj ^ and c- ^ create and annihilate a bosonic 
or fermionic particle with color a on site i, respectively. The 
second term describes an on-site interaction with amplitude 
U± between particles with different colors a ^ (3, with 
^a,i = ^i a^i a^ whcrcas the third term corresponds to the 
on-site interaction between particles with the same color. For 
fermions, the last term vanishes because of the Pauli exclusion 
principle. For femionic alkaline-earth atoms, this description 
is very accurate because, since the interactions between the 
particles do not depend on the nuclear spin, these systems ex- 
hibit a full SU(N) symmetry. In the bosonic case, there is 
in general an additional on-site interaction proportional to the 
square of the total spin that can be antiferromagnetic (as in 
^^Na) or ferromagnetic (as in ^^Rb). This interaction is small 



where S^^ is defined by S^^Ij) = S^^^la). The values of 
the coupling constants Jy and J± depend on the statistics 
of the particles. To second order in t/U\\ and t/U±, they 
are given by Jy = J_\_ = 2t^ /U±_ for fermions, and by 
J|| = 2t^{l/U^ - l/U\\) and J^ = -2t^/U^ for bosons. '^'^^^ 
Thus, the experimentally accessible parameter range is given 
by J|| = J^ > for fermions, and J^ < min(0, — Jy) 
for bosons. Throughout the paper we will use the notation 
J|| = cos and J± = sin to parametrize the Hamiltonian by 
a single parameter 0. 

In the fermionic case with equal couplings Jy = J± 
(6 = 7r/4) the Hamiltonian corresponds to the antiferromag- 
netic (AF) SU(A/^) symmetric Heisenberg model. For N = 2, 
we recover the well-known AF Heisenberg S = 1/2 model. 
The case of larger N has been the subject of many theoreti- 
cal studies, in which a large variety of different ground states 
have been found. For example, for N = 3, 3. three sublattice 
Neel ordered state has been predicted on the square and trian- 
gular lattice, ' whereas on the honeycomb and the kagome 
lattice the ground state can be understood as a generalized 
valence-bond solid state. ^"^^ Even more exotic ground states 
have been predicted for larger N, such as a dimer-Neel or- 
dered state on the square lattice, an algebraic spin-orbital 
liquid on the honeycomb lattice (for A^ = 4), and chiral 
spin liquids. ^"^'^^'^^ 



The anisotropic case of the model of Eq. (2) (Jy ^ J±) 
can be seen as a generahzation of the anisotropic S = 1/2 
Heisenberg (XXZ) model for A/" = 2, for which the phase 
diagram on the square lattice is well known: if | Jy | > \J±\ 
the diagonal terms in the Hamiltonian are dominant and the 
ground state is an antiferromagnetic (ferromagnetic) state for 
J|| > (J|| < 0), whereas if |J||| < \J±\ the off-diagonal 
terms are dominant, and the ground state is a superfluid with 
uniform (2-sublattice) superfluid order parameter for J± < 
(J_\_ > 0). However, for larger N the ground state phase 
diagram is unknown. 

In this article we focus on the case N = 3, and study 
the full ground state phase diagram of the model (2) on a 
square lattice, using linear flavor wave theory (LFWT, Sec. II), 
infinite projected entangled-pair states (iPEPS, Sec. Ill) and 
continuous-time world-line quantum Monte-Carlo (CTQMC, 
Sec. IV). With the latter method we also explore the finite tem- 
perature phase diagram for J^ < 0, in which case there is no 
negative sign problem. 



II. LINEAR FLAVOR WAVE THEORY (LFWT) 
A. The method 

The LFWT is a generalization of the spin-wave expansion 
ofN = 2 models to A^ > 2. The large parameter of LFWT 
is the filling denoted by M (it is usually called spin and de- 
noted by S for N = 2). The representation used on each site 
corresponds to a Young tableau with M boxes arranged hori- 
zontally. 

First, one chooses a classical state \(j)) to which quantum 
fluctuations will be added. These states (the product states, 
described in Sec. II B) are the ground states of the M = oo 
model and are the subject of the next subsection. We then per- 
form a local SU(A/') transformation such that the image state 
10') has M bosons of color a = N on each site. Next the 
transformed Hamiltonian (expressed in terms of new opera- 
tors) is expanded in powers of Xj^fM and its truncation to 
quadratic operators is solved by a Bogoliubov transformation. 
The results for our model are given in Sec. II C. The possible 
phase transitions are discussed in Sec. II D. 

B. Product states 

Let us define aj^ as an operator creating a bosonic parti- 
cle of color Oi Sit site i and |0) as the vacuum of particles. 
A state |(/)i) at site i can be fully described by a unitary N 
component complex vector z^ = (^ii, ^i2, • • • , ^iAr) such that 

\^i) = Y^a ^icx0^ia\^)i- The local U(l) phasc of the z^ vec- 
tor has no importance. Thus the state corresponds to a single 
point of the complex projective space CP(2), but the vectorial 
notation will be conserved to keep notations simple. 

Product states are simply tensor products of single site 
states, i.e. non-entangled states. They are the generalization 
of the classical spin states to the N color models (for TV = 2, 
product states are fully magnetized states characterized by a 



unitary three-dimensional real vector at each site). For a lat- 
tice of Ns sites, a product state is described by A^s(A^ — 1) 
complex numbers, to be compared with (A^^^ — 1) for any 
state of the full Hilbert space. In this reduced space, the 
minimal energy is variational. Since the square lattice with 
nearest-neighbor interactions is bipartite, the lowest energy 
state minimizes the energy of each link {ij) which depends 
only on z^ and Zj : 



^ij = (^11 ~ ^^) 



■ ^^ I X] ^*a 



(3) 



Depending on 0, there are four phases of product ground 
states for A/" = 3 (see the inner circle of Fig. 1): 

• the 2-sublattice superfluid (2subSF) phase (7r/4 < 6 < 
37r/4). There are three distinct ground state mani- 
folds. The states of one of them are characterized by a 
phase (j) such that the vectors on the two sublattices are 

(t^'^T?'^) ^""^ (-TI'^T?'^)- The two remain- 
ing manifolds are obtained by a cyclic permutation of 
the vector components. This leads to a set of ground 
states homeomorph to C3 x U{1). The ground state 
energy per bond is £^2subSF = (^|| - ^±)/2. 

• the antiferromagnetic (AF) phase (arctan(— 1/2) < 
6 < 7r/4): the vectors are chosen among (1,0,0), 
(0,1,0) and (0,0,1) and set on the lattice in such a way 
that two neighboring sites have different vectors. These 
states are the ground states of the classical antiferro- 
magnetic Potts model, which has an infinite degeneracy, 
and a ground state energy per bond of Eaf = 0. 

• the uniform superfluid (USF) phase (— 37r/4 < 6 < 
arctan(2) — 7r/2): all vectors are identical on the lattice 
and the components fulfill \zia\ = l/V^. The ground 
state manifold is homeomorph to a torus (/7(1) x /7(1)). 
The states can be parametrized by z = ( ^ , ^ , ^ ) 



with a -\- t/j -\- (j) 
J||)/3 per bond. 



0, with an energy £^usf = {2Ji 



• the ferromagnetic (F) phase (— 57r/4 < 6 < — 37r/4): 
all vectors are identical on the lattice and equal to ei- 
ther (1,0,0), (0,1,0) or (0,0,1). These states are 
the ground states of the classical ferromagnetic Potts 
model. In fact, these states are the exact quantum 
ground states in this range of 0. The ground state mani- 
fold is homeomorph to C3, and the energy is Ef — J\\ 
per bond. 

Two values of d do not break the SU(3) symmetry (red 
points in Fig. 1): 

• 6> = 7r/4: This is the AF SU(3) symmetric point. Every 
state with perpendicular neighboring vectors is a ground 
state and there exist an extensive number of such states. 
At this point, the 2 sublattice AF and the 2subSF phases 
of Fig. 1 are linked by a global SU(3) transformation. 




FIG. 1: Phase diagram obtained by linear flavor wave theory (in- 
ner circle) and iPEPS (outer circle). The different phases, 2- 
sublattice superfluid phase (2subSF), 3-sublattice antiferromagnetic 
phase (3subAF), 3-sublattice superfluid phase (3subSF), uniform su- 
perfluid (USF), and ferromagnetic (F), are described in the text. The 
two arrows indicate the direction of the displacement of the bound- 
aries when harmonic fluctuations are taken into account. (J||, J±_) = 
(cos^, sin^). 



• 6 = — 37r/4: This is the F SU(3) symmetric point. The 
ground states are the ferromagnetic states, i.e. those 
with uniform z^. The ground state manifold is home- 
omorph to CP{2). At this point, the F and the USF 
phases of Fig. 1 belong to this manifold, and they are 
linked by a global SU(3) transformation. 

For N = 2, since the square lattice is bipartite, we have 
the same properties for as for —0 and this change of param- 
eters can be compensated by a multiplication of the vectors 
of one sublattice by the matrix iaz (rotation of the spins of tt 
around the z axis in the corresponding spin model). It is this 
property that solves the sign problem of the Heisenberg spin 
model on the square lattice. Then = 37r/4 and — 7r/4 are 
too SU(2) symmetric points for A^ = 2. For A^ = 3 we have a 
reduced SU(2) symmetry at these angles: states with one vec- 
tor component equal to zero everywhere and their transforms 
have the same energy for 6 and —6, thus the subset of such 
ground states of a SU(3) symmetric point are sent to degen- 
erate subsets for opposite 0. It gives a continuum of ground 
states for = 37r/4 including both the 2subSF and F states. 
At ^ = — 7r/4, we do not have such a continuum since the z 
vectors of the ground states have three non zero components. 



C. Quantum fluctuations at the harmonic level 

Harmonic quantum fluctuations lift some of the degenera- 
cies found in the product state analysis of the previous section, 
however, not all of them: when degenerated ground states are 
connected via a symmetry of the Hamiltonian, then the har- 
monic (and higher order) fluctuations cannot lift their degen- 



eracy. This is the case at the SU(3) symmetric point of the F- 
USF boundary (6 = — 37r/4), because the two product states 
can be transformed into each other by a global SU(3) transfor- 
mation. 

The degeneracy is not lifted either at the F-2subSF phase 
boundary (6 = 37r/4), but for a different reason: the F and 
2subSF states both belong to a continuum of states linked by 
a set of transformations homeomorph to SU(2), but that are 
not in the group of the Hamiltonian symmetries (the origin of 
this reduced symmetry was given in the previous section). The 
states with vectors (a, 6, 0) and (—a, 6, 0) on both sublattices 
have the same energy for any complex numbers a and b such 
that |ap + |6p = 1. The ground state manifold is thus ex- 
tended to SU(2) and contains both the F and the 2subSF state. 
Since such states are exact quantum ground states (they have 
the same energy per bond as an isolated bond), the degeneracy 
remains at any order. 

At all other boundaries of the LFWT phase diagram and 
for any value of in the AF phase (since there is an exten- 
sive degeneracy of product- states), the term of order M will 
lift the classical degeneracy. To calculate the first order quan- 
tum correction, we first apply a local SU(A^) transformation to 

the operators S^^ that become S^^^ such that the vector z^ is 
transformed into (0, 0, 1). We then perform the following re- 
placement by bosonic operators in the modified Hamiltonian: 



^fNN 






(4) 



M 






(a and /3 are colors different from N) and solve the Hamilto- 
nian on the square lattice neglecting the 0{^/M) terms. 

Among the extensive number of AF states, we will only 
consider the 2 sublattice order (noted 2subAF) and the 3 sub- 
lattice stripe order (noted 3subAF) depicted in Fig. 2. Each 
state of Sec. II B is stable up to the first order correction in 
the interval of where is is the ground state. Moreover, the 
2-3subAF are both stable until = — 7r/4. For a product 
state, stability means that the derived quadratic Hamiltonian 
has a ground state (the Bogoliubov transformation is possi- 
ble). At the classical boundary between two orders, we can 
compare the corrections of order M of both phases and de- 
termine which one has the lowest energy to this order. Since 
simply taking M = 1 does not give an accurate estimate of 
the energy (one would need further orders), one cannot be 
quantitative regarding the boundary, and we just indicate in 
which direction it shifts due to harmonic fluctuations (arrows 
on Fig. 1). 

The corrections for the 2subAF and 2subSF phases are the 
same at the AF SU(3) symmetric point because of the symme- 
tries, but since 2subAF is degenerate with other product states 
with different symmetry properties, the first order corrections 
of these other states will be different and possibly more favor- 
able. Thus, the boundary between AF and 2subSF phases can 
still be pushed towards larger 6>'s. 

We now give the numerical results and the evolutions for all 
boundaries. At the AF-USF boundary, the correction per site 





FIG. 2: The two states considered in the calculation of the first non 
trivial order of LFWT among the infinite number of AF product 
states: 2subAF and 3subAF. 



is -0.12M for the 3subAF, -0.058M for the 2subAF and 
-0.28M for the USF phase. Thus we expect that the USF 
phase is stable until 6 > arctan(2) — 7r/2. 

At the AF-2subSF boundary, the correction is — 0.51M for 
the 3subAF and -0.22M for the 2subAF and 2subSF. It sug- 
gests that the 3subAF survives for 6 > 7r/4. But the degen- 
erate 3subAF states obtained by global SU(3) transformations 
at the SU(3) symmetric point do not have the same classical 
energy as soon as we are not exactly at the = n/A point. 
For smaller angles, the better states are those described pre- 
viously, with vectors among (1,0,0), (0,1,0) and (0,0,1). 
For larger angles, we can reach E = (Jy — J±)/3 on a 

Imk with vectors among (7^,^,75), i^^J^^j'^) 
and i^.f^J^)^ with j = e2W3 and a + V^ + </> = 0. 
It will thus be this 3 sublattice superfluid stripe order dressed 
by quantum fluctuations that will survive for 6 > 7r/4, noted 
3subSF in Fig. 1 (the ground state manifold is homeomorph 
toZ2^ X [7(1)2). 

All these predictions will be confirmed by the iPEPS and/or 
QMC measurements. 



D. Symmetries and phase transitions 

From the broken symmetries of the product- states of 
Sec. II B, we can have an idea of the finite temperature phase 
transitions. Equivalently, we can look at the topology of the 
ground state manifold. The Mermin- Wagner theorem states 
that a connected ground state manifold cannot give rise to a 
finite order phase transition in two dimensions, except at zero 
temperature. This is the case for the USF phase. 

When the ground state manifold consists of several discon- 
nected regions (when some points cannot be linked by a con- 
tinuous path in the manifold), then there is at least one first or 
second order phase transition. The F and 2subSF phases will 
give rise to such a phase transition with a C3 order parameter. 
The 3subAF and 3subSF phases break a P3 x Z2 symmetry, 
that can be restored through one or more transitions (P3 is the 
group of permutations for three elements). Moreover, it is not 
excluded that 2subAF becomes more stable than 3subAF at fi- 
nite temperature, and that it gives rise to its own transitions. 

Besides the finite order phase transitions, infinite order 
phase transitions associated to topological defects can oc- 
cur when the components of the ground state manifold are not 



simply connected (some paths are not continuously shrinkable 
to a point). This is the case for the USF, 2subSF, and 3subSF 
phases, possessing pairs of defects at low temperature. The 
defects encountered in the 2subSF phase are in Z (the first ho- 
motopy group of U{1) is 7ri(/7(l)) = Z), whereas those of 
the USF and 3subSF phases are in 7ri(/7(l) x /7(1)) = Z x Z. 
Thus, a Kosterlitz-Thouless (KT) transition could occur in 
the 2subSF phase and a KT type phase transition in the USF 
and 3subSF phases. The defects encountered in these two 
phases are the same as the dislocations of triangular two di- 
mensional crystals and the phase transition is similar to its 
melting. 

The case of the USF phase is detailed below. Among the 
three phases with defects, USF is the only one in which de- 
fects do not coexist with a discrete symmetry breaking and for 
which the occurence of a topological phase transition makes 
no doubt. In the two other phases, 2subSF and 3subSF, we 
have both defects and discrete symmetry breaking. It can lead 
to a variety of phase transitions, with the possibility of frac- 
tional vortices, as in the case of the AF XY model on the trian- 
gular lattice. '^~^'^ The larger discrete symmetry and the more 
complicated topological defects for the 3subSF phase require 
further investigations about the nature of the phase transitions, 
which go beyond the scope of the present paper. 

We concentrate here on the USF topological phase transi- 
tion. If such a transition occurs in the 3subSF phase, the same 
description can be easily adapted. We neglect the z compo- 
nent fluctuations and construct a subset of product states in- 
teracting via the same Hamiltonian and possessing the same 
superfluid properties. More precisely, we only consider prod- 
uct states whith Zj = ( ^y=- , ^y=- , ^y=- J . We can fix the gauge 
such that aj + ipj -\-(j)j =0 and parametrize a the state of a site 
state with only two real parameters r]j = aj and fij = '^^ ~^^ . 

Since J_\_ < 0, we note J± = —\J±\ so that the sign of the 
different quantities that will follow is obvious. The energy on 
a bond then reads: 
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where Arj = rjj - 
and A/i gives: 

Hi 



- r]i and A/i = fij — jii. An expansion in At/ 
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Note that 



with the energy of a uniform state £^0 = 
a uniform state is labelled by a point (r/, /i) of the order pa- 
rameter space, whose topological properties will be discussed 
below. Taking the continuum limit of this model on the square 
lattice, we obtain 



H 



dr 2E, 



^o + ^((Vr/)^ + (V/i)^)), (7) 



If we look for the local minima of the energy, we obtain the 
two conditions 



n,, 



v'^T] = 0, vV = 0. 



(8) 



These conditions are verified if rj and /i are constant (global 
minimum of energy, 2Eq per surface unit), but also for con- 
figurations with local singularities called vortices (metastable 
states). This model looks very similar to two uncoupled XY 
models Hij = —\J±\ (cos (At/) + cos(A/i)), that also leads to 
Eq. (7) in the continuum limit, but they are different because 
of the nature of the topological defects. The topological de- 
fects are linked to the set of equivalent points in the covering 
group of the order parameter space: the plane (r/, /i). 

In the XY model, the covering group of the order parame- 
ter space [/(I) is the real line R, and two points are equivalent 
if their distance is a multiple of 27r. A closed path in the or- 
der parameter space corresponds to a path in R between two 
equivalent points, characterized by a relative integer n^, the 
distance between these two points divided by 27r, also called 
the winding number. 

In the two uncoupled XY models, the covering group of 
the order parameter space is the real plane R^ , and equivalent 
points are on a square lattice of spacing 27r. A closed path in 
the order parameter space is now characterized by two relative 
integers {rirj.n^). 

In our model (two coupled XY models), the covering group 
of the order parameter is still a plane, but equivalent points 
(t/, /i) form a triangular lattice of spacing 27r (related to the 
Burger's vectors of a triangular crystal ). The most symmet- 
ric way to characterize the topological properties of a closed 
path in the order parameter space is to use three relative inte- 
gers {wi^W2j ws) defined by 
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linked by i^i + 1^2 + ^3 = and represented on Fig. 3. 

We now want to calculate the energy Ey{wi^W2j ws) of a 
configuration with a vortex of winding numbers {wi,W2, ws) 
at position r = and satisfying Eq. (8). Since the problem 
has a spherical symmetry, r]{r) and /i(r) only depend on the 
modulus r. Thus |V/i(r)| = ^ and |V7/(r)| = ^^^^^^. We 
obtain from Eq. (7) 



3r 



Ey{wi,W2,ws) = — ^(^1 +^2 -^wl)\n-. (12) 

We needed two cut-off length: the lattice spacing a and the 
lattice size L. Six elementary vortices exist, with an en- 
ergy f |J±|lnf: (^1,^2,^3) = (1,-1,0), (0,1,-1), 
(—1, 0, 1) and their antivortices with opposite w's. 

We have here detailed the case of a three color USE phase, 
but similar N color USE phases exist. Then, the number of 
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FIG. 3: Covering groups of the order parameter space of different 
models. All points equivalent to the origin are depicted. The winding 
numbers of the green points are indicated by green lines. 



elementary vortices would be A/" — 1. The transition associ- 
ated to these vortices is numerically characterized through the 
winding number (see Sec. IV). 



III. INFINITE PROJECTED ENTANGLED-PAIR STATES 
A. Method 

An infinite projected-entangled pair state (iPEPS) is a varia- 
tional tensor network ansatz for two-dimensional ground state 
wave functions in the thermodynamic limit,'^^''^^ where the ac- 
curacy can be systematically controlled by the so-called bond 
dimension D. It consists of a unit cell of rank-5 tensors with 
one tensor per lattice site, which is periodically repeated on 
the lattice. The first index of a tensor carries the local Hilbert 
space of the corresponding lattice site, and the four remaining 
indices, called auxiliary bonds, connect the tensor to its four 
nearest neighbors. Details on the iPEPS method can be found 
in Refs. 31-33, and also in Ref. 12, where we used the same 
approach to study the AE SU(3) symmetric point of the cur- 
rent model. 

Eor the experts we note that the optimization of the ten- 
sors is done through imaginary time evolution using the so- 
called simple update, ~ and we verified some simulations 
results also with the full update. ^^ The corner-transfer matrix 
method^^'^ ' is used to contract the tensor network, where 
the error of the approximate contraction can be controlled by 
the boundary dimension x- The error due to the finite x is 
small compared to the symbol sizes. To simulate the state with 
3-sublattice (diagonal) order we used tensors with Z^ symme- 
try ' to improve the efficiency. In the superfluid phases we 
cannot exploit this symmetry since it is spontaneously broken 
in the thermodynamic limit. 

By performing simulations with different (rectangular) unit 
cell sizes in iPEPS we can represent states with different type 



of lattice symmetry breaking. For example, a two-sublattice 
state requires a 2 x 2 unit cell with two different tensors A 
and B arranged in a checkerboard order, whereas a three- 
sublattice state requires a 3 x 3 unit cell with three tensors 
arranged in a stripe order. In practice we perform simula- 
tions with different unit cell sizes and check which ansatz 
yields the lowest variational energy. Finally we note that the 
optimization scheme (imaginary time evolution) requires an 
ansatz with at least two different tensors A and B, and this is 
why the smallest unit cell we consider is 2 x 2 (containing two 
different tensors), even for the uniform phases. 



B. Results 

Our iPEPS results are summarized in Fig. 4: In Fig. 4(a) we 
present the local color densities of the three colors, given by 
pie-charts, on each site in the iPEPS unit cell for the different 
states in the phase diagram. Figure 4(b) shows the ground 
state energy obtained for different values of D in the whole 
parameter range of 0, and in Fig. 4(c) we plot the following 
order parameters. 
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averaged over all sites in the unit cell. A finite value of m^ in- 
dicates that there is long-range color (diagonal) order, whereas 
a finite rriod implies off-diagonal long-range order (a super- 
fluid phase). 

In the following we discuss the results for the individual 
phases and then explain how we accurately determined the 
phase boundaries. 

In agreement with the findings in the previous sections, 
iPEPS predicts a ferromagnetic (F) phase between — 57r/4 < 
< — 37r/4. In this phase each site is occupied by the same 
color, as shown in Fig. 4(a). Since the ferromagnetic state is 
a product state it can be trivially represented by an iPEPS, 
i.e. with a bond dimension D = 1. This is why the energies 
obtained with D = 1 and the D = 8 shown in Fig.4(b) are 
identical in this phase. 

In the range -37r/4 < < -0.109(3)7r iPEPS confirms 
the uniform superfluid phase (USE). The states in this phase 
do not have color order, as can be seen in Fig.4(a) where all 
color densities on all sites are the same. However, we find a 
finite, uniform off-diagonal order parameter rriod, indicating a 
superfluid phase. 

In the range -0.109(3)7r < < 0.48(l)7r we clearly 
find the lowest variational energy with a 3-sublattice iPEPS 
ansatz. This range also includes the SU(3) symmetric point 
= 7r/4, where iPEPS previously found a 3-sublattice Neel 
ordered state. For -0.109(3)7r < 6 < l/47r we find a 3- 
sublattice diagonal order (3subAF in Fig. 4(a) with = 0.27r) 
where each site is occupied by one dominant color, and the 
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FIG. 4: iPEPS results as a function of 0. (a) Color densities in 
the iPEPS unit cell in the different phases, (b) Energy per site for 
different values of D. (c) Plot of the order parameters rrid and 
rriod obtained with D — 8 which indicate long-range color order 
and/or (off-diagonal) superfluid order, respectively, (d) First order 
phase transition between the USF and the SsubAF phase occurring 
for 0/n = —0.109(3). (e) First order phase transition between the 
SsubSF and the 2subSF phase occurring for 0/n ^ 0.48(1). 



three colors are arranged in a stripe pattern. The three sub- 
lattice stripe structure persists for 6 > 7r/4, however, in this 
range it is the offdiagonal terms in rriod which are finite and 
the diagonal terms vanish (3subSF phase). 

We also confirm the 2subSF phase with iPEPS in the range 
0.48(l)7r < 6 < 3/47r. Fig. 4(a) shows that the state has 
a uniform diagonal order (each site is equally populated by 
only two colors with the third color being absent), however, 
the off-diagonal terms exhibit a two-sublattice structure (not 
shown). 

Next we discuss the phase transition between the USF and 
3subAF phase. We find clear evidence of a first order phase 
transition, with a jump of the order parameter at the transition 
value. In order to accurately determine the critical value 6c we 
make use of the hysteresis effect in the vicinity of a first order 
phase transition: When we initialize a state in the USF phase 
and change ^ to a value slightly above 6c then the state will 
remain in the USF phase (because it is metastable). Similarly 
we can obtain the energy of the 3subAF state slightly below 
6c. We then determine 6c = —0.109(3) from the intersection 
of the energies of the two states, as shown in Fig.4(d). The er- 
ror bar of the phase boundary includes the uncertainties when 
extrapolating the energies to the infinite D limit. 

In a similar way we have determined the location of the 
first order transition 6c = 0.48 (l)7r between the 3subSF and 
the 2subSF phase, as shown in Fig. 4(e). 

We did not find that the phase boundaries of the ferromag- 
netic phase change when going to larger bond dimensions. 
Thus, the location of the phase transitions is the same as pre- 
dicted classically, i.e. 6 = — 37r/4 and 6 = 37r/4. 



IV. CONTINUOUS TIME WORLD LINE MONTE-CARLO 
WITH CLUSTER UPDATES 

The iPEPS numerical simulations give access to infor- 
mations at zero temperature. Using quantum Monte Carlo 
(QMC) methods, finite temperature observables can be mea- 
sured. Moreover, QMC methods work in the full Hilbert space 
and thus is a non biased numerical method. They give inter- 
esting complementary informations to the iPEPS and LFWT 
results. 

Unfortunately, the so-called sign problem occurs as soon as 
some non-diagonal elements of the Hamiltonian are positive. 
To study the model of Eq. (2), we chose a basis of states con- 
sisting of the product states <S)i\ai) = |Qfo, • • • , o^n^-i), where 
ai is the color of the boson on site i and Ns is the number of 
sites of the lattice. These are eigenstates of all the /Sf ". The 
sign of the non-diagonal elements of the Hamiltonian matrix 
is the sign of Jj_ . The sign problem only allows us to simulate 
this model for J^ < 0. 

In some specific cases, the sign problem can be overcome. 
On an infinite or open one dimensional chain, the physical 
properties of the N color model are unchanged by the trans- 
formation 6^—6, whatever N. This allows us to use the 
QMC method whatever the value of 6, although mainly the 
SU(A^) symmetric point was considered. ' In two dimen- 
sions, the same transformation still conserves the physical 
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TABLE I: Fabrication rules for clusters. First column: all possible 
graphs G that can be placed on a bond of the lattice at an imaginary 
time. Three following columns: Ag(c), which is equal to 1 when 
G is accepted and otherwise. A dash (— ) means that the graph is 
incompatible with the local configuration c. Last two columns: time 
constant used to determine when a graph is inserted on a bond with 
a Poisson law. e is a small number ensuring the ergodicity (e was set 
to 0.01 in our simulations). 



properties for A^ = 2 (see discussion in Sec. II C), but not 
any more f or A^ > 2, where qualitatively different phases can 
be obtained for two opposite values of 6, as can be seen in the 
phase diagram of Fig. 1. QMC is thus limited to J^ < 0. 



A. The algorithm 

We use a continuous time world-line algorithm with cluster 
updates,"^^'"^"^ adapted to our 3-color model. We express the 
partition function Z as a path integral over the configurations 
(j) : r -^ ^(t), where r is the imaginary time going from 
to /3 = 1/T (T is the temperature) and (j){r) is a basis state. 
The functions (j) that contribute to Z can be represented by 
^(0) and by a set of world line crossings {(z, j, r)} at which 
the colors of two neighboring sites i and j at a time r are 
exchanged. A local configuration c on a link ij at time r is 
represented by 



_^a,(r+)a,(r+) 
ai{r )aj{r ) 



(15) 



Cluster algorithms are well known for 2-color (spin) models. 
Here we generalize them by chosing randomly two different 
colors p and q out of the 3 and by constructing clusters on 
which only these two colors are encountered. The steps to 
construct a cluster are the following. We first randomly place 
the graphs depicted in the first column of Tab. I on each link 
of the configuration using a Poisson distribution whose time 
constant W{G) is given in the last columns of Tab. I and ac- 
cepting them if Ag(c) = 1 (if a color being neither p nor q 
appears in the local configuration, the graph is rejected). Then 
we assign graphs to the world-line crossings between p and q 
colors with a probability proportional to W{G). At the places 
where no graph has been attributed, we follow the path with 
the same color. Finally we follow each constructed cluster and 
exchange p and q on it with a probability 1/2. This completes 
one Monte Carlo step. During the whole simulation, n steps 
are performed. 



We use periodic square lattices of linear size L. The goal is 
now to calculate averages of observables such as the correla- 
tions over the n sampled configurations. 

We define the diagonal correlation on site j as 
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with Yj the position vector of site j. This structure factor is 
normalized in such a way that ^ ^j^ C'(k) = 1. 

B. Results 

We have studied the sector —it/ 2 < < Oof Fig. 1 us- 
ing the previously detailed algorithm. Confirming LFWT and 
iPEPS results, we found two different phases at low temper- 
ature: the USF and the 3subAF stripe order. We studied the 
thermal phase transitions generated by these two phases (see 
Fig. 5). For USF, it is the KT type phase transition linked to 
the winding number W described in Sec. II D. The 3subAF 
stripe order gives rise to a first order transition due to the 
Ps X Z2 symmetry breaking. 

The clearest sign of the P3 x Z2 symmetry breaking is 
the structure factor in the reciprocal space. Bragg peaks ap- 
pear at two points of coordinates q = ±(27r/3, 27r/3) or 
±(27r/3, — 27r/3) depending on the direction of the stripes 
(see Fig. 6(c)). Their evolution with the size of the lattice as 
I/^ indicates long range order. The transition temperature in- 
creases with more negative O's (see Fig. 5) and corresponds to 
a first order transition: coexistence of ordered and unordered 
phases occurs during the simulations. When 6 :^ — O.IItt, no 
trace of the 3subAF phase is found at any temperature. At 
large temperature, the structure factors have peaked maxima 
at point q = (tt, tt) (similar to Fig. 6(b)), but the finite size 
effects indicate no long range order associated to these wave 
vectors. 

For the KT type phase transition of the USF phase, we 
looked at the evolution of the stiffness T = {W'^)/2T with 
the temperature. The critical temperature is deduced from 
the intersection of the curves with the line T/tt, ' extrap- 
olated to the thermodynamic limit. At = — 7r/2 it is 
Tkt = 0.564(2) (see Fig. 7). Even if this phase cannot be 
characterized by the structure factor as there is no long range 
order, it looks quite different from the high temperature phase. 
The peak at q = (tt, tt) is much broader, as illustrated on 
Fig. 6(a). 
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FIG. 5: Phase diagram for -7r/2 <0 <0 obtained from QMC. The 
dashed line is a KT phase transition, whereas the continuous line is a 
first order transition. 
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FIG. 6: Structure factors obtained on a L = 12 square lattice for the 
model of Eq. (2) with A^ = 3. The white square is the Brillouin zone. 



presented the full phase diagram as a function of the coupling 
parameter 0, where Jy = cosO and Jj_ = sinO are the ex- 
change amplitudes between respectively identical and differ- 
ent colors. The model corresponds to the infinite repulsion 
limit of the fermionic Hubbard model in the = 7r/4 case, 
and of the bosonic Hubbard model for — tt < < — 7r/4. 
Three complementary methods have been used that give 
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V. RESULTS AND CONCLUSIONS 

We have studied a quantum version of the N = 3-state 
Potts model on the square lattice which can be seen as a gen- 
eralization of the well known 5 = 1/2 XXZ model. We have 



FIG. 7: Stiffness T versus temperature for different linear lattice 
sizes L for 6 — — 7r/2 (precision of 0.001: smaller than the symbol 
size). The temperature of crossing with the dashed red line T = T/tt 
is shown in the inset with the interpolation giving a critical tempera- 
ture of T^t = 0.564(2). 
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FIG. 8: Stiffness T versus temperature for different and a lattice 
linear size of L = 12. The temperature of crossing with the dashed 
red line T = T/tt is shown in the inset and is compatible with a 
critical temperature reaching zero for ^ 0.1 Itt. When not visible, 
the error bars are smaller than the symbol size. 



coherent results: linear flavor wave theory, iPEPS and quan- 
tum Monte Carlo. The resulting ground state phase diagram 
as a function of 6 is given in Fig. 1 and includes 5 phases: 
three superfluid phases (USF, 2subSF, and 3subSF), a ferro- 
magnetic phase (F), and a three-sublattice color-ordered phase 
(3subAF). The temperature dependent properties have been 
determined for— 7r/2 < 6 < and are presented in the phase 
diagram of Fig. 5. The 3subAF phase gives rise to a first or- 



der phase transition and the USF phase to a KT type phase 
transition. 

The phase diagram for A^ = 3 has many similarities with 
the one from the S = 1/2 XXZ model (N = 2): A uni- 
form superfluid and a A/^-sublattice superfluid characterized 
by A^ — 1 phase angles, a ferromagnetic phase that can take N 
different colors, and a A'-sublattice antiferromagnetic, stripe- 
ordered phase. Each of these phases for A" = 3 can be seen as 
a generalization of the corresponding phase for A^ = 2. How- 
ever, the A^ = 3 case has one additional phase: the 2subSF 
phase, which is not strictly speaking a generalization of a 
phase of the A" = 2 case, but is in fact equivalent to the 2- 
sublattice superfluid phase of the A^ = 2 case (since only two 
colors are present in this phase). Thus it is interesting that for 
A^ = 3 a 2-color superfluid is competing with a 3 -color su- 
perfluid, and that both phases are stable in a certain parameter 
range. One finds also a competition between a 2-sublattice 
and a 3-sublattice color ordered state (which are degenerate at 
the classical level). In this case, however, it is the 3-sublattice 
color-ordered state which has always a lower energy than the 
2-sublattice state. 

Finally, coming back to possible experimental realizations, 
three phases of this phase diagram should be accessible with 
cold atoms loaded in optical lattices: the USF phase and the 
F phase with bosonic particles, and the 3subAF phase at the 
symmetric AF point with fermionic particles. 
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